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I. SOMVBY 


A study has been completed of mathematically proper boundary conditions for 
unique numerical solution of internal, viscous, subsonic flows in the SSME. The 
study has concentrated on well-posed considerations, with emphasis on 
canputationa 1 efficiency and numerically stable boundary condition statements. 
The method of implementing the established boundary conditions is applicable to 
a wide variety of finite difference and finite element codes, as demonstrated. 
The results of this study are reported herein. 

II. TECHNICAL DISCUSSION 

A. Introduction 

Over the past several years a focus at NASA Marshall Space Flight Center has 
been adaptation and application of canputationa 1 fluid dynamics (CFD) analysis 
techniques to flowfield prediction in components of the SSME. Several 
"olympiads" have been held, wherein purveyors of CFD codes have developed and 
compared solutions for model problem definition analyses to the turn-around 
duct-transfer duct SSME geometry. The SSME geometry is defined to these codes 
via construction of meshes that possess boundary segments roughly coincident 
with solid walls and containing convenient flew inlet and outlet planes. The 
numerical simulation of the associated flowfield is defined via appropriate 
specification of constraints on the (Navier-Stokes) conservation law system 
variables, e.g., velocity and pressure, over the entire boundary of the mesh. 

This study examined boundary condition specifications for the CFD models, 
with emphasis on mathematical well-posedness with physical consistency. The 
SSME flowfield is characterized as complex turbulent three-dimensional and 
unsteady, at high Reynolds number and lew subsonic Mach number, ie., essentially 
incompressible. Mathematically, the CFD algorithms/codes applied to this 
problem definition fall into two distinguishable categories. One family (GIM, 
PAGE, VAST, ...} sets up solution as the unsteady evolution of a hyperbolic 
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conservation law system with the assumption of a compressible fluid satisfying a 
polytropic gas law statement. Conversely, the second family (INS3D, PHOENIX, 
SIMPLE, FIDAP) specifically assumes an incompressible fluid, and directly seeks 
the steady-state solution without specifying a (physically significant) equation 
of state. The PHOENIX and SIMPLE algorithms seek the steady-state through a 
pressure relaxation procedure that explicitly requires pressure boundary 
condition specifications. FIDAP uses a finite element penalty method to totally 
replace the appearance of pressure. Alternatively, the INS3D theory employs a 
pseudo-compressibility concept, yielding a hyperbolic conservation law-appearing 
statement for pressure that (only) requires approximation of the normal pressure 
gradient at boundaries. Thus, the mathematical boundary condition aspect of 
INS3D is more analogous to that of GIM, et.al., than PHOENIX, et.al. 

The following sections examine CFD algorithm boundary condition issues from 
the standpoint of, a) hyperbolic conservation law systems, and, b) pressure 
relaxation in an incompressible flew specification. 

B. Hyperbolic Conservation law Algorithms 

The conservation law system governing the kinematics, kinetics and 
thermodynamics of a viscous, heat-conducting fluid is generally termed the 
Navier-Stokes equations. The Cartesian tensor indicial form is, 

t, . ap a 

<p * * * *. < v ) = 0 

a(pu.) d 

= ~r + * <v“, + - v = ° 
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rt \ _ a(pe) a 
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where i> is density, Pir^ is the momentum vector, p is pressure, 5— is the 

Kronecker delta, e is specific total energy, and y is the ratio of specific 

heats for a polytropic gas law fluid. The expression of constitutive properties 

of the fluid is contained in the stress tensor and heat flux vector qj. For 

simple fluids and laminar flew, the accepted forms are, 

du. du n du, 1 

— + — - — - — - 6 .. 

to. dx. 3 dx k uj ( 5 ) 

8 T 

( 6 ) 

l 

where the dynamic (molecular) viscosity P(T), and the thermal oonductibity k(T), 
are weak functions of temperature T, and the Reynolds number is Re=(pUL/u) ref. 

In the limit of large Reynolds number, an inviscid, non-heat conducting 
assumption renders equations 5-6 identically zero. The resulting form of eqn 1- 
4 is termed the Euler equations, a homogeneous hyperbolic conservation law 
system. Alternatively, enforcing a statistical averaging procedure yields a 
Reynolds-averaged Navier-Stokes system that introduces the concept of a Reynolds 
stress tensor and additional governing partial differential equations, eg., the 
two-equation turbulent kinetic energy-isotropic dissipation function system. 

In either instance the generic form of the governing equation system is, 

a, < 

-+-♦« = <> (7) 

J 

where q contains the Navier-Stokes/Euler/Reynolds-averaged Navier-Stokes 
dependent variable set, f j is the corresponding flux vector and s is a 
source/sink term, eg.. 
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The single-point closure equations for c^j, qj, kj and e j are (cf.. Baker, 
1986), “ 


o. = o. — pa a 
U u * J 


q = q + p H'u' — u' a — u a. 

v j i u ‘ u 

* = (c. piFiF - - ls )± 

j \ k * j e Re y / as. 

I 

e = (c pu' u' . - ) — - 
' V * K * J e / a*. 


(9) 


The Euler equation form is contained in eqn 22 by deletion of k, e , kj, e j, qj, 
j and s, and replacement of pH byp e + yp* 

A familiar alternative form for eqn 7 is established following imposition of 
a coordinate transformation x^ = x^(rij)» where j = ££, n, C ) is any 
(curvilinear) coordinate system. One particularly useful form is to align the 
coordinate 5 with the direction of principal flow, whereupon eqn 7 can be 
written as 
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where E, F and G are the (Euler) flux vector components, E v , F v and G v are the 
constitutive closure model components (containing o^j, qj, kj and e j). Both are 
expressed in terms of scalar components in the n j coordinate system, ie.. 
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where J is the Jacobian of the coordinate transformation, and the convection 
velocity (0) vector contravariant scalar components parallel to the (K , n ) 
coordinate system are, 

U = t, u + t, v + S, u> 

x y z 

V = nu + n u + n u> 

x l y l t 

w = (,u + + C. f w ( 12 ) 


The constitutive scalar components E v , F v and G v are formed in the similar 
manner. 

The conservation law system, eqn 7 or 10, is either hyperbolic or an 
initial-valued, incompletely elliptic boundary value problem, dependent upon the 
constitutive closure definitions (inviscid, visoous/turbulent). The solution 
domain ft x t is a bounded subregion of a portion of the SSME duct region, and 
the boundary conditions on 3ft , mathematically consistent with a well-posed 
problem, have been examined by Dutt(1985) following a dependent variable 
transformation to "entropy variables." The entropy transformation of the 
primitive (Euler) dependent variable set q = (p , p u^, pe), eqn 1-3, has been 
examined, cf., Harten(1981), Osher, et.al. (1984). For a family of strictly 
convex entropy functions, Mai let (1985) and others show that for the (Euler) 
extension to viscous and heat-conducting fluids, the sole "useful" entropy 
function is the thermodynamic entropy p s. Hence, the transformation is, 

V(q) = -ps = -p log (p/f!> 6 ) (13) 

and the entropy flux functions are 

f j = - pu^s) = - m i (os) (14) 

The transformation to V(q) symmetrizes the conservation law statement, eqn 1-3, 
yielding a nonlinear energy estimate (for sufficiently smooth solutions to the 
mixed initial-boundary value problem) that corresponds to the Clausius-Duhem 
inequality (second law of thermodynamics). 
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For the (Euler equation)' hyperbolic conservation law form, the appropriate 
number of boundary conditions on oo is (Strickwerda, 1977): supersonic inflow 
(5), subsonic inflow (4), supersonic outflow (none), subsonic outflow (1). 
For the Navier Stokes equations, 5 (4) boundary conditions are required at 
inflow (outflow). Dutt(1985) develops the set of "maximal dissipative" boundary 
conditions for eqn 1-6 for the (Navier-Stokes) boundary condition statement 
form, 

£R || + Sq - b (15) 

where R is a matrix of rank at most 4, 5 is the coordinate normal to'3!), and for 
the Euler definition ( e = 0), Sq=b is a proper form for the (unperturbed Euler) 
hyperbolic problem. The derived outflow boundary conditions are inner (dot) 
products of eqn 5-6 with the outward pointing unit normal vector hj, ie., 

- a^ = b ± , 1 < i < 3 

q-j.fij = 0 (16) 

where is the velocity contravariant scalar component parallel to £. Any SSME 
application involves only subsonic outflow, hence a- L (= aC ) > 0 and b^ is a 
constant and all other components of a^ and b^ vanish as does the normal heat 
flux. The derived inflow boundary condition couples the influx definition pU-^ = 
b Q with the more general form of eqn 16. 

PUi = b 0 

ai j .n j - a i U^ = bj^ , 1 _< i < 3 (17) 

- q-j.nj " a 4 T = b 4 

In eqn 17, the subscript bar on the contravariant velocity vector denotes the 
index is not summed, and a^ and b^ (0 < i < 4) are constants subject to 
constraints. For the SSME definition of subsonic inflow (M < 1), a^ = 0 = b-^ 
while a 2 , > b Q /2 and a 4 > 1/2 ( - y + 2 ) b q + e, where e corresponds to 

the inverse Reynolds number, see eqn 5. Of sane interest, the last expression 
in eqn 17 is directly amenable to physical interpretation; in expanded form, 
using eqn 6, 
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(18) 


-qj.fij = k-~ .fij = KVT.n = b 4 + a 4 T = h (T-T r ) 

Thus, a 4 is interpreted as the boundary heat convection coefficient h, and T r is 
the heat exchange reference temperature. 

Equations 15-18, in concert with Strickwerda's constraints, encompass the 
range of boundary conditions, that are mathematics 1 ly well-posed and numerically 
stable, for SSME flew CFD simulations formulated as approximate procedures for 
solving the hyperbolic conservation law Euler extension to Navier-Stokes. Since 
SSME flcwfields are uniformly subsonic, then only one exit Dirichlet boundary 
condition is allowed, taken as the static pressure. Up to four inlet Dirichlet 
boundary conditions are mathematically permitted; however, eqn 17 suggests that 
all but the mass influx be replaced by Neumann constraints. Further, if the 
mass influx is specified, then the inlet pressure may not be specified (unless a 
region of supersonic flow exists between inlet and outlet). Conversely, a 
Dirichlet specification of inlet total pressure could be made, whereupon the 
mass flux will become determined by the f lowpath total pressure drop. The 
reported SSME simulations using PAGE and VAST have generally employed the 
former, while the INS3D simulations have generally used the latter. 

C. Incompressible Navier-Stokes Algorithms 

As noted at the end of Section A, the alternative SSME simulation CFD 
construction class utilizes incompressibility directly to recast the Navier- 
Stokes system into a nan-hyperbolic conservation law form. In these procedures, 
eg., PHOENIX, SIMPLE and FIDAP, the pressure distribution is derived indirectly 
from the velocity constraint of divergence-freeness, hence no equation of state 
is (need be) assumed to exist. Thus, no (Dirichlet) pressure boundary 
conditions are needed or appropriate in defining the CFD simulation, although it 
is well published that PHOENIX and SIMPLE employ a "pressure correction 
equation" to achieve convergence to a numerical approximation of divergence- 
freeness. This is not strictly exact, as will be developed. 
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The mentioned class of incompressible Navier-Stokes algorithms can be viewed 
in a unified manner as decisions made in evaluating a Taylor series on the time 
evolution of the velocity field u fe , where boldface defines a vector field. 
Assuming knowledge of the solution at time t^, where t n+1 = + At, we have, 

uP +1 = up + Atu£ + ... (19) 

The incompressible form of eqn 2, plus eqn 5, provides the expression for the 
time-derivative in eqn 19, hence, 

u n+ ^ = t/ 1 - At [u n «vtf 1 + V p — Re - -*- V ] + ... (20) 

The basic CFD algorithm theoretical choice lies in selection of vp in eqn 
20. A finite element penalty algorithm replaces the variable p with the 
approximation to continuity, 

p = - XV .u (21) 

where X is a large 0(10^ ) constant. Ihe superscript tilde denotes that u is an 
approximation to a divergence-free velocity field. Finally, eqn 21 is evaluated 
at the new time t n+ ^, hence eqn 20-21 is an implicit expression. 

Hie basic theory for the SIMPLE-class of incompressible Navier-Stokes CFD 
algorithms is similarly developed directly from eqn 20. If ^Vp n is employed, 
then the solution u n+ ^ does not satisfy the continuity equation. Hence, define 
a new pressure p n+1 that produces (assumption) a divergenoe-free velocity field 
u n+1 . Writing eqn 20 for both pressures, taking the sum and subtracting yields, 

Vx (if 1+1 - uP +1 ) = 0 (22) 

Thus, the distinction between the predicted and the continuity-satisfying 

velocity fields at t n+ ^ can at most be the gradient of a scalar field Q, ie., 

^n+l _ jjp+l = (23) 

Subtracting eqn 22 into the incompressible form of eqn 1 yields 

V 2 d> = - V ,u n+1 (24) 

The boundary condition for eqn 23, for the harmonic function 0, is obtained frcm 
eqn 22 as 

V<t>.n = (u 0 * 1 - u" +1 ) .n (25) 
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where n is the unit outwarding pointing vector normal to the solution domain o. 
Once dP is determined, using an appropriate (CFD) algorithm for eqn 23-24, then 
the "corrected pressure" field is 

P" +1 - P" - (26) 

At steady-state convergence, eqn 23 becomes homogeneous, hence |d) h l - O for a 
(single) Dirichlet boundary specification in concert with eqn 24. Thus in the 
limit, p is the pressure field that coexists with a computed approximation to 
the divergence-free velocity field u* 1 . 

Viewing eqn 20, 22-25, there is no admissible pressure boundary condition 
specification for the SIMPLE-class, CFD incompressible Navier-Stokes algo- 
rithms. Equations 16-17 remain appropriate for the remaining variables, and the 
mass flowrate (eqn 17a) must be specified to create a SSME problem statement. 
The auxiliary variable d) carries the retaining boundary condition specification, 
and eqn 24 is homogeneous Neumann everywhere that the distinction between u* 1 and 
u n must vanish, eg., inlet, solid (no-slip) walls, symmetry planes, etc. The 
sole Dirichlet constraint for (D* 1 therefore can only be applied at a location on 
the mass efflux boundary segment of oo. Assuming the CFD iteration sequence is 
convergent, eqn 25 yields the corresponding static pressure distribution to 
within an arbitrary constant, which can be specified (for example) to match p to 
an inlet or an outlet pressure level. 
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D. Concluding Discussion 

The literature contains numerous results documenting the robustness of the 
CFD algorithm boundary conditions developed in the preceeding sections. The 
compressible hyperbolic conservation law formulation is exhaustively examined in 
Dutt(1985); Chang and Kwak(1984) document the total pressure specification 
option for the INS3D algorithm. To our knowleddge, the interpretation of 
SIMPLE-type algorithms developed in Section C is not cannon knowledge. The 0* 1 
construction for incompressible parabolic Navier-Stokes algorithms is well 
established; Baker(1983, Ch. 6-7) fully documents the range of application of 
the solution statement given as eqn 23-24. 
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